Thư viện

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(stringr)
library(stats)
library(boot)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## 
## The following object is masked from 'package:boot':
## 
##     melanoma
## 
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(leaps)
library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loaded glmnet 4.1-8
library(Metrics)
## 
## Attaching package: 'Metrics'
## 
## The following objects are masked from 'package:caret':
## 
##     precision, recall
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:boot':
## 
##     logit
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(visreg)
library(corrplot)
## corrplot 0.94 loaded
library(lmPerm)
library(dplyr)
library(mgcv)
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.

Đọc và xử lí dữ liệu

data <- read.csv("cleaned_fifa_eda_stats.csv")
head(data)
##       id              name age nationality overall potential
## 1 158023          L. Messi  31   Argentina      94        94
## 2  20801 Cristiano Ronaldo  33    Portugal      94        94
## 3 190871         Neymar Jr  26      Brazil      92        93
## 4 193080            De Gea  27       Spain      91        93
## 5 192985      K. De Bruyne  27     Belgium      91        92
## 6 183277         E. Hazard  27     Belgium      91        91
##                  club value  wage preferred_foot international_reputation
## 1        FC Barcelona 110.5 0.565           Left                        5
## 2            Juventus  77.0 0.405          Right                        5
## 3 Paris Saint-Germain 118.5 0.290          Right                        5
## 4   Manchester United  72.0 0.260          Right                        4
## 5     Manchester City 102.0 0.355          Right                        4
## 6             Chelsea  93.0 0.340          Right                        4
##   weak_foot skill_moves attacking_work_rate deffensive_work_rate body_type
## 1         4           4              Medium               Medium      Lean
## 2         4           5                High                  Low    Stocky
## 3         5           5                High               Medium      Lean
## 4         3           1              Medium               Medium      Lean
## 5         5           4                High                 High    Normal
## 6         4           4                High               Medium    Normal
##     position jersey_number     joined loaned_from contract_valid_until height
## 1    Forward            10 2004-07-01 Not on loan           2021-12-31     67
## 2    Forward             7 2018-07-10 Not on loan           2022-12-31     74
## 3    Forward            10 2017-08-03 Not on loan           2022-12-31     69
## 4 Goalkeeper             1 2011-07-01 Not on loan           2020-12-31     76
## 5 Midfielder             7 2015-08-30 Not on loan           2023-12-31     71
## 6    Forward            10 2012-07-01 Not on loan           2020-12-31     68
##   weight crossing finishing heading_accuracy short_passing volleys dribbling
## 1    159       84        95               70            90      86        97
## 2    183       84        94               89            81      87        88
## 3    150       79        87               62            84      84        96
## 4    168       17        13               21            50      13        18
## 5    154       93        82               55            92      82        86
## 6    163       81        84               61            89      80        95
##   curve fk_accuracy long_passing ball_control acceleration sprint_speed agility
## 1    93          94           87           96           91           86      91
## 2    81          76           77           94           89           91      87
## 3    88          87           78           95           94           90      96
## 4    21          19           51           42           57           58      60
## 5    85          83           91           91           78           76      79
## 6    83          79           83           94           94           88      95
##   reactions balance shot_power jumping stamina strength long_shots aggression
## 1        95      95         85      68      72       59         94         48
## 2        96      70         95      95      88       79         93         63
## 3        94      84         80      61      81       49         82         56
## 4        90      43         31      67      43       64         12         38
## 5        91      77         91      63      90       75         91         76
## 6        90      94         82      56      83       66         80         54
##   interceptions positioning vision penalties composure marking standing_tackle
## 1            22          94     94        75        96      33              28
## 2            29          95     82        85        95      28              31
## 3            36          89     87        81        94      27              24
## 4            30          12     68        40        68      15              21
## 5            61          87     94        79        88      68              58
## 6            41          87     89        86        91      34              27
##   sliding_tackle gk_diving gk_handling gk_kicking gk_positioning gk_reflexes
## 1             26         6          11         15             14           8
## 2             23         7          11         15             14          11
## 3             33         9           9         15             15          11
## 4             13        90          85         87             88          94
## 5             51        15          13          5             10          13
## 6             22        11          12          6              8           8
##   release_clause
## 1          226.5
## 2          127.1
## 3          228.1
## 4          138.6
## 5          196.4
## 6          172.1

Ta loại bỏ các cột không đóng góp giá trị như id, name, nationality, loan_from, club

- id, name: do là các giá trị riêng biệt

- club : do các biến này có thể đóng góp trong thực tế (ví dụ như cầu thủ của Real Madrid thì tiềm năng hơn) tuy nhiên ta sẽ loại bỏ nó vì có quá nhiều câu lạc bộ sẽ làm tăng quá nhiều độ phức tạp của mô hình, tương tự với ‘nationality’ và ‘loan_from’.

data <- data|> dplyr::select(-c(id,name))

data <- data |> dplyr::select(-c(loaned_from,club,nationality))

Xử lí dữ liệu:

- Ta chỉnh các biến định tính về dạng ‘factor’

- Các biến có dạng mm-dd-yyyy thì ta chỉ giữ lại năm (yyyy).

data_numeric <- data |> dplyr::select(where(is.numeric))

all_colname <- names(data)

category_col <- setdiff(all_colname, names(data_numeric))

for (i in 1:length(category_col)) {
  data <- data |>
    mutate(across(all_of(category_col[i]), as.factor))
}

data$contract_valid_until <- as.numeric(substr(data$contract_valid_until, 1, 4))
data$joined <- as.numeric(substr(data$joined, 1, 4))

Tách dữ liệu làm hai tập train test

set.seed(100)
id_train <- sample(1:nrow(data),size = as.integer(0.8* nrow(data)))
id_test <- setdiff(1:nrow(data), id_train)
data_train <- data[id_train,]
data_test <- data[id_test,]

I. Linear Regression cho Overall

1. Mô hình hồi quy với tất cả các biến

model <- lm(formula = overall ~ ., data = data_train)
result_display <- function(model, data_test, target){
  
  # Phần dư và các chỉ số trên tập huấn luyện
  residuals_train <- residuals(model)
  mae_train <- mean(abs(residuals_train))
  mse_train <- mean(residuals_train^2)
  adjr2_train <- summary(model)$adj.r.squared
  
  # Dự đoán trên tập kiểm tra
  y_predict <- predict(model, newdata = data_test)
  
  # Tính các chỉ số trên tập kiểm tra
  mse_test <- mean((data_test[[target]] - y_predict)^2)
  mae_test <- mean(abs(data_test[[target]] - y_predict))

  SSE <- sum((data_test[[target]] - y_predict)^2)
  SST <- sum((data_test[[target]] - mean(data_test[[target]]))^2)  
  
  r_squared_test <- 1 - (SSE / SST)

  n_test <- length(data_test[[target]])  # số mẫu trong tập kiểm tra
  k <- length(model$coefficients) - 1    # số tham số trong mô hình (không tính hằng số)
  
  adj_r_squared_test <- 1 - ((1 - r_squared_test) * (n_test - 1) / (n_test - k - 1))
  
  # Tạo bảng kết quả
  results <- data.frame(
    Metric = c("MSE", "MAE", "Adjusted R2"),
    Train = c(mse_train, mae_train, adjr2_train),
    Test = c(mse_test, mae_test, adj_r_squared_test)
  )
  
  return(results)
}
result_display(model,data_test,'overall')
##        Metric    Train      Test
## 1         MSE 3.105704 3.1407562
## 2         MAE 1.376685 1.3953978
## 3 Adjusted R2 0.936822 0.9334985

2. Lựa chọn mô hình

2.1 Step_wise với kFolds Cross-Validaiton

library(leaps)
set.seed(100)
k = 5 
folds<-sample(rep(1:k,length=nrow(data_train)))
 predict.regsubsets <- function(object, newdata, id_model){
   form <- as.formula(object$call[[2]])
   x_mat <- model.matrix(form, newdata)
   coef_est <- coef(object, id = id_model)
   x_vars <- names(coef_est)
   res <- x_mat[, x_vars] %*% coef_est
   return(as.numeric(res))
 }

Ta sẽ dùng phương pháp là ‘backward’ để giảm thiểu thời gian cho việc chọn biến vì nếu như dùng ‘exhaustive’ cho quá nhiều biến thì rất tốn thời gian và thậm chí không tìm ra.

cv_error_fifa_rj<-matrix(0,nrow=k,ncol=dim(data_train)[2])
for (r in 1: k){
  train <- data_train[folds !=r,]
  validation <- data_train[folds ==r,]
  out_subset_fifa_folds<-regsubsets(x=overall~ .,data=train,method="backward",nvmax=dim(train)[2])
   for(j in 1:dim(train)[2]){
     pred_rj<-predict(out_subset_fifa_folds, newdata=validation,id_model=j)
     cv_error_fifa_rj[r,j]<-(mean((validation$overall-pred_rj)^2))
 }
}
 cv_error_fifa<-colMeans(cv_error_fifa_rj)
ggplot(data= data.frame(x= c(1:dim(data_train)[2]),y=cv_error_fifa),aes(x=x,y=y)) +
 geom_point() +
 geom_line()+
 labs(x="Number of predictors",y="RMSE") +
 theme_bw()

cv_error_fifa
##  [1] 27.593196  7.339429  5.373930  5.149759  4.855027  4.404690  4.130525
##  [8]  3.936668  3.884836  3.804003  3.717270  3.631348  3.550559  3.468331
## [15]  3.412661  3.349174  3.327005  3.318421  3.295094  3.278514  3.254145
## [22]  3.240575  3.229941  3.214501  3.208693  3.203589  3.198037  3.188483
## [29]  3.184845  3.174307  3.165740  3.158478  3.155494  3.149582  3.147407
## [36]  3.147207  3.145309  3.145252  3.143827  3.141378  3.139023  3.138328
## [43]  3.136577  3.136931  3.137426  3.137658  3.138497  3.139476  3.139936
## [50]  3.139705  3.139229  3.139325  3.139597

Dưới đây là mô hình tốt nhất được lựa chọn từ phương pháp này

out_subset_fifa_2 <- regsubsets(x = overall ~ ., data = data_train,method = "backward",nvmax = dim(train)[2])
a <- names(coef(out_subset_fifa_2, id = which.min(cv_error_fifa)))
data_train_subset <- data_train|> dplyr::select(-c(wage,preferred_foot,weak_foot,attacking_work_rate,deffensive_work_rate,body_type,position,joined,height,dribbling,curve,fk_accuracy,long_passing,jumping,long_shots,interceptions))

Ta thấy các thông số đo lường tuy có giảm nhưng không quá đáng kể, nhưng ta đã giảm được độ phức tạp của thuật toán khá nhiều khi bỏ đi một số các biến đầu vào.

model <- lm(formula = overall ~ ., data = data_train_subset)
result_display(model,data_test,'overall')
##        Metric     Train      Test
## 1         MSE 3.2069351 3.2131129
## 2         MAE 1.3995615 1.4078213
## 3 Adjusted R2 0.9348659 0.9324005

2.2 Ridge và Lasso

x <- model.matrix(overall ~ ., data = data_train)[,-1]
y <- data_train$overall
library(glmnet)
out_cv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, type.measure = "mse", nfolds = 5,family = "gaussian")
lambda_grid <- 10^seq(from = 10, to =-2, length = 100)

beta_lambda_lasso <- out_cv_lasso$lambda.min
out_lasso_md <- glmnet(x = x, y = y, alpha = 1,lambda = lambda_grid, family = "gaussian")
predict(out_lasso_md, s = beta_lambda_lasso, type = "coefficients")
## 58 x 1 sparse Matrix of class "dgCMatrix"
##                                        s1
## (Intercept)                 -1.528883e+02
## age                          4.889756e-01
## potential                    4.681959e-01
## value                        1.179531e-01
## wage                         .           
## preferred_footRight          .           
## international_reputation    -3.668166e-01
## weak_foot                    5.964669e-03
## skill_moves                  7.541796e-01
## attacking_work_rateLow       3.901313e-01
## attacking_work_rateMedium    .           
## deffensive_work_rate Low     3.059105e-01
## deffensive_work_rate Medium -3.518333e-02
## body_typeNormal              6.995978e-02
## body_typeStocky              9.750410e-02
## positionForward             -1.001122e+00
## positionGoalkeeper           .           
## positionMidfielder          -9.541762e-01
## jersey_number               -9.447815e-03
## joined                       5.150668e-04
## contract_valid_until         7.033457e-02
## height                       .           
## weight                       8.433614e-03
## crossing                     2.241263e-03
## finishing                    1.150634e-02
## heading_accuracy             3.790878e-02
## short_passing                4.478217e-02
## volleys                      .           
## dribbling                    4.584999e-05
## curve                        .           
## fk_accuracy                  .           
## long_passing                 1.388231e-03
## ball_control                 6.799523e-02
## acceleration                 1.411227e-02
## sprint_speed                 2.390074e-02
## agility                      4.878108e-03
## reactions                    1.347731e-01
## balance                     -1.506928e-02
## shot_power                   1.439248e-02
## jumping                      2.016065e-03
## stamina                      2.921972e-02
## strength                     2.308028e-02
## long_shots                   .           
## aggression                   1.696301e-03
## interceptions                .           
## positioning                 -1.391864e-02
## vision                      -5.738963e-03
## penalties                   -3.289357e-03
## composure                    3.749191e-02
## marking                      6.715413e-03
## standing_tackle              4.937880e-04
## sliding_tackle              -6.323383e-05
## gk_diving                    4.587257e-02
## gk_handling                  3.642575e-02
## gk_kicking                   7.611226e-03
## gk_positioning               2.387573e-02
## gk_reflexes                  5.079425e-02
## release_clause               .

Cross-Validation cho mô hình đã bỏ đi các biến theo Lasso

data_train_lasso <- data_train|>dplyr::select(-c(release_clause,interceptions,long_shots,curve,fk_accuracy,volleys,weight,attacking_work_rate,wage,preferred_foot))

mse_list <- c()

# kfolds = 5
for (i in 1: 5){
  train <- data_train_lasso[folds != i,]
  validate <- data_train_lasso[folds == i,]
  
  fifa_lasso <- lm(formula = overall ~.,data = train)
  
  mse <- (mean((validate$overall-predict(newdata = validate,fifa_lasso))^2))
  mse_list <- c(mse_list,mse)
}
mse_cv <- mean(mse_list)
cat("RMSE: ",mse_cv)
## RMSE:  3.160872

Các thông số bằng phương pháp này cho kết quả tốt hơn so với phương pháp trước nên vì thế ta sẽ dùng mô hình này cho các bước tiếp theo.

data_train_lasso <- data_train|>dplyr::select(-c(release_clause,interceptions,long_shots,curve,fk_accuracy,volleys,weight,attacking_work_rate,wage,preferred_foot))

model <-lm(formula = overall ~.,data = data_train_lasso)

result_display(model,data_test,'overall')
##        Metric     Train     Test
## 1         MSE 3.1328745 3.178213
## 2         MAE 1.3833992 1.402746
## 3 Adjusted R2 0.9363221 0.932931

3. Chuẩn đoán mô hình

3.1 Kiểm tra tính tuyến tính của mô hình

model <- lm(data = data_train_lasso, formula = overall ~. )
library(ggplot2)

ggplot(model,aes(x = .fitted, y = .resid))+
  geom_point()+
  geom_smooth(method = 'loess',se = F)+
  geom_hline(yintercept = 0,linetype = 'dashed')+
  theme_bw()+
  labs(x = "Fitted values", y = "Residuals")
## `geom_smooth()` using formula = 'y ~ x'

Ta thấy dựa vào hình thì giá trị thặng dư khá tốt ở lúc đầu tuy nhiên phần cuối lại có dấu hiệu là đường cong nên vì thế ta sẽ hiệu chỉnh sau.

3.2 Kiểm tra tính tuyến tính từng phần

terms_md <- predict(model, type = "terms")
partial_resid_md <- residuals(model, type = "partial")

col = names(data.frame(terms_md))

# Duyệt qua tất cả các cột trong data và vẽ biểu đồ cho từng cột
for (i in 1:length(col)){

  # Tạo dataframe cho mỗi cột
  data_part_resid_md <- tibble(
  data_ = data_train[, col[i]],
  terms = terms_md[, col[i]],
  partial_resid = partial_resid_md[, col[i]]
  )
  
  # Vẽ biểu đồ và hiển thị biểu đồ cho từng cột
  p <- ggplot(data_part_resid_md, mapping = aes(data_, partial_resid)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
  geom_line(aes(x = data_, y = terms), color = "blue") +
  labs(x = col[i], y = "Partial Residuals") +
  theme_bw()

# In ra biểu đồ
  print(p)
}
## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0004
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 5.9127e-15
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0804

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2020
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

3.3 Kiểm tra đồng nhất phương sai

ggplot(model, aes(.fitted, sqrt(abs(.stdresid)))) +
  geom_point(na.rm = TRUE) +
  geom_smooth(method = "loess", na.rm = TRUE, se = FALSE) +
  labs(x = "Fitted Values", y = expression(sqrt("|Standardized residuals|"))) +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

Từ hình ảnh ta thấy các điểm đầu khá chuẩn chỉ có những điểm dữ lệu cuối có xu hướng đường cong rõ rệt tương tự như phần 3.1

3.4 Kiểm tra điểm ngoại lai

ggplot(model, aes(.hat, .stdresid)) +
 geom_point(aes(size = .cooksd)) +
 xlab("Leverage") + ylab("Standardized Residuals") +
 scale_size_continuous("Cook's Distance", range = c(1, 6)) +
 theme_bw() +
 theme(legend.position = "bottom")

std_resid_md <- rstandard(model)
hat_values_md <- hatvalues(model)
cooks_D_md <- cooks.distance(model)
 data_cooks_md <- tibble(id_point = 1:nrow(data_train_lasso),
 rstand = std_resid_md, hats = hat_values_md,
 cooks = cooks_D_md, overall = data_train_lasso$overall)
 data_cooks_md <- data_cooks_md |> arrange(desc(cooks))
 head(data_cooks_md)
## # A tibble: 6 × 5
##   id_point rstand   hats   cooks overall
##      <int>  <dbl>  <dbl>   <dbl>   <int>
## 1      879  -4.38 0.0399 0.0170       94
## 2     6110  -4.47 0.0336 0.0148       91
## 3     8803  -5.28 0.0239 0.0145       91
## 4     4357  -4.89 0.0195 0.0101       94
## 5    13229  -3.52 0.0224 0.00606      91
## 6     3209  -4.81 0.0106 0.00528      56

Ta sẽ thực hiện loại bỏ Outlier cho mô hình

Những điểm nào có cook distance lớn hơn 4/n với n là tổng số quan sát thì ta sẽ loại bỏ những điểm này

id_exclude <- (data_cooks_md |> filter(cooks > 4/nrow(data_train_lasso)))$id_point
id_rest <- setdiff(data_cooks_md$id_point,id_exclude)
data_train_lasso <- data_train_lasso[id_rest,]
model <- lm(data = data_train_lasso, formula = overall ~. )


result_display(model,data_test,'overall')
##        Metric     Train      Test
## 1         MSE 2.2461048 3.2155173
## 2         MAE 1.2121397 1.3954416
## 3 Adjusted R2 0.9514614 0.9321437

Ta thấy mô hình cho kết quả tốt hơn ở tập train so với ban đầu. Vì các điểm Outlier đã bị loại bỏ giúp cho mô hình học tốt hơn. Tuy nhiên ở tập test đã tăng lên. Ta sẽ giải quyết vấn đề này ở phần 3

3.5 Kiểm tra đa cộng tuyến

model <- lm(data = data_train_lasso, formula = overall ~. )
library(car)
df_vif <- data.frame(vif(model))|>arrange(desc(GVIF)) |>filter(GVIF >5)

df_vif
##                        GVIF Df GVIF..1..2.Df..
## position         340.117025  3        2.642032
## gk_diving         29.814994  1        5.460311
## gk_reflexes       29.737673  1        5.453226
## gk_handling       27.429321  1        5.237301
## gk_positioning    26.026394  1        5.101607
## standing_tackle   24.929892  1        4.992984
## gk_kicking        23.005440  1        4.796399
## sliding_tackle    22.969645  1        4.792666
## ball_control      19.398961  1        4.404425
## dribbling         16.015987  1        4.001998
## short_passing     13.112069  1        3.621059
## finishing          9.762605  1        3.124517
## positioning        9.463701  1        3.076313
## acceleration       8.902934  1        2.983779
## long_passing       7.514671  1        2.741290
## sprint_speed       7.368550  1        2.714507
## marking            6.802705  1        2.608200
## heading_accuracy   6.732468  1        2.594700
## crossing           6.111110  1        2.472066
## shot_power         5.123720  1        2.263564
remove_multicollinearity <- function(feature, data_train) {
  return (data_train |> dplyr::select(-all_of(feature)))
}
while(TRUE){
  model <- lm(data = data_train_lasso, formula = overall ~. )
  
  df_vif <- data.frame(vif(model))|>arrange(desc(GVIF)) |>filter(GVIF >5)|>slice_max(GVIF, n = 1)
  vec_feature <- row.names(df_vif)
  
  
  if (length(vec_feature) == 0){
    break
  }
  
  data_train_lasso <- remove_multicollinearity(vec_feature,data_train_lasso)
}
model <- lm(data = data_train_lasso, formula = overall ~. )


result_display(model,data_test,'overall')
##        Metric     Train     Test
## 1         MSE 3.0568351 4.060881
## 2         MAE 1.3910795 1.550162
## 3 Adjusted R2 0.9340255 0.914720

Sau khi loại bỏ hết đa cộng tuyến thì mô hình trở nên rất tê, tuy nhiên nó đã giảm được độ phức tạp cuẩ thuật toán khá nhiều

4. Mở rộng mô hình

upgrade_model <- lm(formula = overall ~ poly(age,4)+poly(value,2) +poly(potential,2)+., data = data_train_lasso)

result_display(upgrade_model,data_test,'overall')
##        Metric     Train      Test
## 1         MSE 1.5574229 2.6839806
## 2         MAE 0.9411552 1.0527663
## 3 Adjusted R2 0.9663734 0.9434984

Ta sẽ sử dụng hàm mũ cho các biến như ‘age’, ‘value’, ‘potential’ ta thấy mô hình đã cải thiện rõ rệt và cho thông số rất tốt.

data_tst <- data_test|>dplyr::select(names(data_train_lasso))
new_obs <- data_tst[10,]
new_obs
##    age overall potential value international_reputation weak_foot skill_moves
## 57  26      86        87    52                        3         4           4
##    deffensive_work_rate body_type jersey_number joined contract_valid_until
## 57               Medium      Lean            10   2016                 2023
##    height crossing heading_accuracy long_passing sprint_speed agility reactions
## 57     69       73               62           71           93      91        86
##    balance shot_power jumping stamina strength aggression vision penalties
## 57      86         82      75      84       67         73     82        71
##    composure marking
## 57        80      42

5. Xây dựng khoảng tin cậy

5.1 Xay dựng khoảng tin cậy cho từng biến

fun_boot_md <- function(data, ind, formula, ...){
 data_new <- data[ind,]
 out_md <- lm(formula = formula, data = data_new, ...)
 return(out_md$coefficients)
}
set.seed(100)

out_boot_md <- boot(data = data_train_lasso, statistic = fun_boot_md, R = 1000,
 formula = overall ~ poly(age,4)+ poly(value,2)+ poly(potential,2)+ .,parallel = "multicore")
pvals <- sapply(1:ncol(out_boot_md$t),function(x) {
 qt0 <- mean(out_boot_md$t[, x] <= 0)
 if (!is.na(qt0)){
   if (qt0 < 0.5) {
   return(2*qt0)
   } else {
   return(2*(1- qt0))
   }
 }
 }
 )
pvals <- unlist(pvals)


cat(pvals, sep = " ")
## 0.674 0 0 0 0 0 0 0 0 0 0.386 0 0.102 0.538 0 0 0 0.23 0 0.002 0 0 0 0 0 0 0 0.072 0 0.748 0 0.164 0.112 0 0 0

Bên trên và p_value cho từng biến

5.2 Khoảng tin cậy cho trung bình

fun_model_predict_mu <- function(data,ind,formula,new_observation){
  data_new <- data[ind,]
  model <- lm(formula = formula,data = data_new)
  pred <- predict(model,newdata = new_observation)
  return (pred)
}
library(boot)
set.seed(100)

out_boot_pred <- boot(data = data_train_lasso, statistic = fun_model_predict_mu,formula = overall ~ poly(age,4)+ poly(value,2)+ poly(potential,2)+.,new_observation = new_obs, R = 1000, parallel = "multicore")
df <- data.frame(t = out_boot_pred$t)

# Create the histogram using ggplot
ggplot(df, aes(x = t)) +
  geom_histogram(binwidth = 0.2, fill = "lightblue", color = "black") +
  labs(x = "Bootstrap Estimates", y = "Frequency") +
  theme_minimal()

quantile(out_boot_pred$t,c(0.025,0.975))
##     2.5%    97.5% 
## 82.71068 85.24096

5.3 Khoảng tin cậy cho biến dự đoán

 resid <- residuals(upgrade_model)
 y_pd_pci <- out_boot_pred$t + sample(resid, size = 1000, replace = TRUE)
 quantile(y_pd_pci, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 81.41045 86.92959

II. Linear Regression cho Potential

1. Mô hình hồi quy với tất cả các biến

model <- lm(formula = potential ~ ., data = data_train)
result_display(model,data_test,'potential')
##        Metric     Train      Test
## 1         MSE 5.9263871 5.8796982
## 2         MAE 1.9126669 1.9026382
## 3 Adjusted R2 0.8443281 0.8349253

2. Lựa chọn mô hình

2.1 Step_wise và kFolds Cross-Validaiton

library(leaps)
set.seed(101)
k = 5 
folds<-sample(rep(1:k,length=nrow(data_train)))
 predict.regsubsets <- function(object, newdata, id_model){
   form <- as.formula(object$call[[2]])
   x_mat <- model.matrix(form, newdata)
   coef_est <- coef(object, id = id_model)
   x_vars <- names(coef_est)
   res <- x_mat[, x_vars] %*% coef_est
   return(as.numeric(res))
 }
cv_error_fifa_rj<-matrix(0,nrow=k,ncol=dim(data_train)[2])
for (r in 1: k){
  train <- data_train[folds !=r,]
  validation <- data_train[folds ==r,]
  out_subset_fifa_folds<-regsubsets(x=potential~ .,data=train,method="backward",nvmax=dim(train)[2])
   for(j in 1:dim(train)[2]){
     pred_rj<-predict(out_subset_fifa_folds, newdata=validation,id_model=j)
     cv_error_fifa_rj[r,j]<-(mean((validation$potential-pred_rj)^2))
 }
}
 cv_error_fifa<-colMeans(cv_error_fifa_rj)
ggplot(data= data.frame(x= c(1:dim(data_train)[2]),y=cv_error_fifa),aes(x=x,y=y)) +
 geom_point() +
 geom_line()+
 labs(x="Number of predictors",y="RMSE") +
 theme_bw()

out_subset_fifa_2 <- regsubsets(x = potential ~ ., data = data_train,method = "backward",nvmax = dim(data_train)[2])
a <- names(coef(out_subset_fifa_2, id = which.min(cv_error_fifa)))
setdiff(names(data_train),a)
##  [1] "potential"            "wage"                 "preferred_foot"      
##  [4] "skill_moves"          "attacking_work_rate"  "deffensive_work_rate"
##  [7] "body_type"            "position"             "weight"              
## [10] "heading_accuracy"     "short_passing"        "volleys"             
## [13] "dribbling"            "curve"                "fk_accuracy"         
## [16] "ball_control"         "reactions"            "shot_power"          
## [19] "jumping"              "aggression"           "interceptions"       
## [22] "gk_diving"            "gk_kicking"
data_train_subset <- data_train|> dplyr::select(-c(
  wage, preferred_foot, skill_moves, attacking_work_rate, 
  deffensive_work_rate, body_type, position, weight, heading_accuracy, 
  dribbling, curve, fk_accuracy, reactions, shot_power, jumping, 
  interceptions, gk_kicking)
)
model <- lm(formula = potential ~ ., data = data_train_subset)
result_display(model,data_test,'potential')
##        Metric     Train      Test
## 1         MSE 5.9605624 5.9274720
## 2         MAE 1.9178459 1.9110744
## 3 Adjusted R2 0.8436898 0.8346958

2.2 Ridge và Lasso

x <- model.matrix(potential ~ ., data = data_train)[,-1]
y <- data_train$potential
library(glmnet)
out_cv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, type.measure = "mse", nfolds = 5,family = "gaussian")
lambda_grid <- 10^seq(from = 10, to =-2, length = 100)

beta_lambda_lasso <- out_cv_lasso$lambda.min
out_lasso_md <- glmnet(x = x, y = y, alpha = 1,lambda = lambda_grid, family = "gaussian")
predict(out_lasso_md, s = beta_lambda_lasso, type = "coefficients")
## 58 x 1 sparse Matrix of class "dgCMatrix"
##                                        s1
## (Intercept)                 -3.638345e+01
## age                         -9.240860e-01
## overall                      8.707526e-01
## value                       -1.681563e-02
## wage                         .           
## preferred_footRight         -3.750818e-02
## international_reputation     1.002604e+00
## weak_foot                    3.983963e-02
## skill_moves                  .           
## attacking_work_rateLow       1.120708e-02
## attacking_work_rateMedium   -8.281785e-04
## deffensive_work_rate Low    -1.386661e-01
## deffensive_work_rate Medium -6.593011e-02
## body_typeNormal             -2.084916e-01
## body_typeStocky             -2.174540e-01
## positionForward              3.466959e-01
## positionGoalkeeper          -5.170995e-01
## positionMidfielder          -1.499088e-02
## jersey_number                1.200514e-02
## joined                      -3.078260e-02
## contract_valid_until         6.910695e-02
## height                      -3.453569e-02
## weight                      -2.620536e-03
## crossing                    -1.238857e-02
## finishing                   -1.430527e-04
## heading_accuracy             .           
## short_passing                .           
## volleys                     -1.254810e-03
## dribbling                    .           
## curve                        .           
## fk_accuracy                 -2.195305e-04
## long_passing                -2.264711e-03
## ball_control                 .           
## acceleration                -4.019492e-03
## sprint_speed                -1.835308e-02
## agility                     -3.280621e-03
## reactions                    .           
## balance                      7.270086e-03
## shot_power                   .           
## jumping                      .           
## stamina                     -4.183314e-02
## strength                    -1.330754e-02
## long_shots                  -7.347920e-03
## aggression                   .           
## interceptions                .           
## positioning                 -6.392837e-03
## vision                       1.304703e-02
## penalties                    1.525798e-02
## composure                    2.325342e-02
## marking                      8.231183e-03
## standing_tackle              .           
## sliding_tackle               7.824327e-03
## gk_diving                    .           
## gk_handling                  .           
## gk_kicking                   .           
## gk_positioning               3.407980e-03
## gk_reflexes                 -3.800857e-03
## release_clause               .

Cross- Validation cho Lasso Model

data_train_lasso <- data_train|>dplyr::select(-c(release_clause,gk_reflexes,gk_kicking,gk_handling,gk_diving,sliding_tackle,interceptions,shot_power,reactions,ball_control,long_passing,fk_accuracy,curve,dribbling,short_passing,finishing,skill_moves,wage))

mse_list <- c()

# kfolds = 5
for (i in 1: 5){
  train <- data_train_lasso[folds != i,]
  validate <- data_train_lasso[folds == i,]
  
  fifa_lasso <- lm(formula = potential ~.,data = train)
  
  mse <- (mean((validate$potential-predict(newdata = validate,fifa_lasso))^2))
  mse_list <- c(mse_list,mse)
}
mse_cv <- mean(mse_list)
cat("RMSE: ",mse_cv)
## RMSE:  6.015211
model <-lm(formula = potential ~.,data = data_train_lasso)

result_display(model,data_test,'potential')
##        Metric     Train      Test
## 1         MSE 5.9694848 5.9365298
## 2         MAE 1.9179766 1.9108106
## 3 Adjusted R2 0.8434086 0.8342419

Ta chọn bộ data ‘data_train_subset’ cho những phần tiếp theo.

3. Chuẩn đoán mô hình

3.1 Kiểm tra tính tuyến tính của mô hình

model <- lm(data = data_train_subset, formula = potential ~. )
library(ggplot2)

ggplot(model,aes(x = .fitted, y = .resid))+
  geom_point()+
  geom_smooth(method = 'loess',se = F)+
  geom_hline(yintercept = 0,linetype = 'dashed')+
  theme_bw()+
  labs(x = "Fitted values", y = "Residuals")
## `geom_smooth()` using formula = 'y ~ x'

Tính tuyến tính của mô hình là không thỏa chủ yếu ở những điểm đầu và cuối.

3.2 Kiểm tra tính tuyến tính từng phần

terms_md <- predict(model, type = "terms")
partial_resid_md <- residuals(model, type = "partial")

col = names(data.frame(terms_md))

# Duyệt qua tất cả các cột trong data và vẽ biểu đồ cho từng cột
for (i in 1:length(col)){

  # Tạo dataframe cho mỗi cột
  data_part_resid_md <- tibble(
  data_ = data_train[, col[i]],
  terms = terms_md[, col[i]],
  partial_resid = partial_resid_md[, col[i]]
  )
  
  # Vẽ biểu đồ và hiển thị biểu đồ cho từng cột
  p <- ggplot(data_part_resid_md, mapping = aes(data_, partial_resid)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
  geom_line(aes(x = data_, y = terms), color = "blue") +
  labs(x = col[i], y = "Partial Residuals") +
  theme_bw()

# In ra biểu đồ
  print(p)
}
## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0004
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 5.9127e-15
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0804

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0

## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2020
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

3.3 Kiểm tra đồng nhất phương sai

ggplot(model, aes(.fitted, sqrt(abs(.stdresid)))) +
  geom_point(na.rm = TRUE) +
  geom_smooth(method = "loess", na.rm = TRUE, se = FALSE) +
  labs(x = "Fitted Values", y = expression(sqrt("|Standardized residuals|"))) +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

Mô hình không đồng nhất phương sai

3.4 Kiểm tra điểm ngoại lai

ggplot(model, aes(.hat, .stdresid)) +
 geom_point(aes(size = .cooksd)) +
 xlab("Leverage") + ylab("Standardized Residuals") +
 scale_size_continuous("Cook's Distance", range = c(1, 6)) +
 theme_bw() +
 theme(legend.position = "bottom")

std_resid_md <- rstandard(model)
hat_values_md <- hatvalues(model)
cooks_D_md <- cooks.distance(model)

Ta sẽ thực hiện loại bỏ Outlier cho mô hình

Những điểm nào có cook distance lớn hơn 4/n với n là tổng số quan sát thì ta sẽ loại bỏ những điểm này

 data_cooks_md <- tibble(id_point = 1:nrow(data_train_subset),
 rstand = std_resid_md, hats = hat_values_md,
 cooks = cooks_D_md)
 data_cooks_md <- data_cooks_md |> arrange(desc(cooks))
 head(data_cooks_md)
## # A tibble: 6 × 4
##   id_point rstand    hats   cooks
##      <int>  <dbl>   <dbl>   <dbl>
## 1     4357   3.41 0.0434  0.0147 
## 2     8052   2.22 0.0503  0.00724
## 3    12025   1.79 0.0461  0.00429
## 4     3485   3.47 0.0105  0.00355
## 5     2317   5.17 0.00284 0.00211
## 6     1716   4.49 0.00368 0.00206
id_exclude <- (data_cooks_md |> filter(cooks > 4/nrow(data_train_subset)))$id_point
id_rest <- setdiff(data_cooks_md$id_point,id_exclude)
data_train_subset <- data_train_subset[id_rest,]
model <- lm(data = data_train_subset, formula = potential ~. )

result_display(model,data_test,'potential')
##        Metric     Train      Test
## 1         MSE 4.3800636 6.0144317
## 2         MAE 1.6926823 1.9039860
## 3 Adjusted R2 0.8786617 0.8322707

Ta thấy mô hình cho kết quả tốt hơn ở tập train so với ban đầu. Vì các điểm Outlier đã bị loại bỏ giúp cho mô hình học tốt hơn. Tuy nhiên ở tập test đã tăng lên. Ta sẽ giải quyết vấn đề này ở phần 3

3.5 Kiểm tra đa cộng tuyến

model <- lm(data = data_train_subset, formula = potential ~. )
library(car)
data.frame(vif = vif(model))|>arrange(desc(vif))
##                                 vif
## value                    101.815814
## release_clause            95.415415
## gk_reflexes               27.058978
## gk_diving                 26.881912
## gk_positioning            24.627297
## gk_handling               24.619208
## standing_tackle           24.408502
## sliding_tackle            22.429390
## ball_control              16.145992
## short_passing             13.199222
## finishing                 10.640687
## positioning                8.979948
## acceleration               8.723754
## long_passing               7.297106
## sprint_speed               7.285258
## long_shots                 7.202759
## marking                    6.633039
## overall                    6.346339
## volleys                    6.256924
## crossing                   5.177654
## agility                    4.672618
## balance                    4.498720
## penalties                  4.434403
## vision                     4.170814
## stamina                    3.830155
## aggression                 3.630437
## composure                  3.620232
## height                     3.391360
## strength                   2.776957
## international_reputation   2.088940
## age                        2.011132
## weak_foot                  1.178769
## contract_valid_until       1.173539
## jersey_number              1.108738
## joined                     1.103443
while(TRUE){
  model <- lm(data = data_train_subset, formula = potential ~. )
  
  df_vif <- data.frame(vif = vif(model))|>arrange(desc(vif))|>filter(vif >5)|>slice_max(vif, n = 1)
  vec_feature <- row.names(df_vif)
  
 
  if (length(vec_feature) == 0){
    break
  }
 
  data_train_subset<- remove_multicollinearity(vec_feature,data_train_subset)
}
model <- lm(data = data_train_subset, formula = potential ~. )
result_display(model,data_test,'potential')
##        Metric    Train      Test
## 1         MSE 4.450271 6.1205107
## 2         MAE 1.705565 1.9174774
## 3 Adjusted R2 0.876844 0.8299836

4. Mở rộng mô hình

upgrade_model <- lm(formula = potential ~ poly(age,4)+., data = data_train_subset)

result_display(upgrade_model,data_test,'potential')
##        Metric     Train      Test
## 1         MSE 2.1285713 3.2472175
## 2         MAE 1.0426589 1.1800429
## 3 Adjusted R2 0.9410803 0.9096891
summary <- summary(upgrade_model)
cat("Adjusted R2 train:",summary$adj.r.squared)
## Adjusted R2 train: 0.9410803
data_tst <- data_test|>dplyr::select(names(data_train_subset))
new_obs <- data_tst[4,]
new_obs
##    age overall potential international_reputation weak_foot jersey_number
## 39  33      88        88                        3         2             1
##    joined contract_valid_until height crossing volleys long_passing
## 39   2012                 2021     76       12      12           34
##    sprint_speed agility balance stamina strength aggression vision penalties
## 39           55      47      36      41       71         25     41        23
##    composure marking release_clause
## 39        69      25             51

5. Xây dựng khoảng tin cậy

5.1 Xay dựng khoảng tin cậy cho từng biến

fun_boot_md <- function(data, ind, formula, ...){
 data_new <- data[ind,]
 out_md <- lm(formula = formula, data = data_new, ...)
 return(out_md$coefficients)
}
set.seed(100)

out_boot_md <- boot(data = data_train_subset, statistic = fun_boot_md, R = 1000,
 formula = potential ~ poly(age,4)+ .,parallel = "multicore")
pvals <- sapply(1:ncol(out_boot_md$t),function(x) {
 qt0 <- mean(out_boot_md$t[, x] <= 0)
 if (!is.na(qt0)){
   if (qt0 < 0.5) {
   return(2*qt0)
   } else {
   return(2*(1- qt0))
   }
 }
 }
 )
pvals <- unlist(pvals)


cat(pvals, sep = " ")
## 0.09 0 0 0 0 0 0 0.508 0 0 0 0.898 0 0.006 0.982 0.014 0.022 0.028 0 0.2 0 0.47 0.002 0.35 0 0.368

5.2 Khoảng tin cậy cho trung bình

fun_model_predict_mu <- function(data,ind,formula,new_observation,knot){
  data_new <- data[ind,]
  model <- lm(formula = formula,data = data_new)
  pred <- predict(model,newdata = new_observation)
  return (pred)
}
library(boot)
set.seed(100)

out_boot_pred <- boot(data = data_train_subset, statistic = fun_model_predict_mu,formula =
                       potential ~ poly(age,4)+.,new_observation = new_obs, R = 1000, parallel = "multicore")
df <- data.frame(t = out_boot_pred$t)

# Create the histogram using ggplot
ggplot(df, aes(x = t)) +
  geom_histogram(binwidth = 0.02, fill = "lightblue", color = "black") +
  
  labs(x = "Bootstrap Estimates", y = "Frequency") +
  theme_minimal()

quantile(out_boot_pred$t,c(0.025,0.975))
##     2.5%    97.5% 
## 88.19318 88.50122

5.3 Khoảng tin cậy cho biến dự đoán

 resid <- residuals(upgrade_model)
 y_pd_pci <- out_boot_pred$t + sample(resid, size = 1000, replace = TRUE)
 quantile(y_pd_pci, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 85.18643 92.16681

III. Linear Regresion cho Value

simplified_data <- data|> dplyr::mutate(wage = wage * 1000)
value_lm <- lm(value ~ ., data = simplified_data)

summary(value_lm)
## 
## Call:
## lm(formula = value ~ ., data = simplified_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5958 -0.1200 -0.0012  0.0994 12.2424 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  6.212e+00  8.374e+00   0.742 0.458256    
## age                         -1.910e-02  2.186e-03  -8.737  < 2e-16 ***
## overall                      3.463e-02  2.492e-03  13.894  < 2e-16 ***
## potential                   -1.887e-02  1.812e-03 -10.415  < 2e-16 ***
## wage                         6.375e-03  4.121e-04  15.471  < 2e-16 ***
## preferred_footRight         -2.878e-03  1.115e-02  -0.258 0.796337    
## international_reputation     3.996e-01  1.621e-02  24.657  < 2e-16 ***
## weak_foot                    1.405e-02  7.369e-03   1.906 0.056669 .  
## skill_moves                  7.028e-03  1.206e-02   0.583 0.560214    
## attacking_work_rateLow      -2.762e-04  2.425e-02  -0.011 0.990913    
## attacking_work_rateMedium   -8.047e-03  1.164e-02  -0.691 0.489405    
## deffensive_work_rate Low    -2.003e-03  2.107e-02  -0.095 0.924250    
## deffensive_work_rate Medium  1.367e-03  1.309e-02   0.104 0.916853    
## body_typeNormal              1.765e-02  1.038e-02   1.700 0.089120 .  
## body_typeStocky              3.999e-02  2.122e-02   1.884 0.059518 .  
## positionForward              2.842e-02  2.551e-02   1.114 0.265218    
## positionGoalkeeper           4.991e-01  9.608e-02   5.194 2.08e-07 ***
## positionMidfielder           1.527e-02  1.795e-02   0.851 0.394799    
## jersey_number               -5.115e-04  2.921e-04  -1.751 0.079930 .  
## joined                       1.660e-03  2.223e-03   0.746 0.455377    
## contract_valid_until        -4.942e-03  3.711e-03  -1.332 0.183021    
## height                      -4.716e-03  3.468e-03  -1.360 0.173931    
## weight                       7.156e-04  5.250e-04   1.363 0.172832    
## crossing                    -2.233e-03  6.336e-04  -3.525 0.000425 ***
## finishing                   -3.162e-03  7.689e-04  -4.112 3.94e-05 ***
## heading_accuracy            -9.633e-04  6.773e-04  -1.422 0.154999    
## short_passing               -2.754e-03  1.085e-03  -2.539 0.011135 *  
## volleys                      3.783e-03  6.564e-04   5.763 8.39e-09 ***
## dribbling                   -6.404e-04  9.505e-04  -0.674 0.500492    
## curve                       -1.662e-03  6.483e-04  -2.564 0.010364 *  
## fk_accuracy                  1.958e-03  5.756e-04   3.402 0.000672 ***
## long_passing                 2.529e-03  7.985e-04   3.168 0.001540 ** 
## ball_control                -1.806e-03  1.184e-03  -1.525 0.127287    
## acceleration                 6.035e-04  8.892e-04   0.679 0.497316    
## sprint_speed                -9.413e-05  8.250e-04  -0.114 0.909165    
## agility                      7.978e-04  6.603e-04   1.208 0.226968    
## reactions                    1.772e-03  1.006e-03   1.761 0.078335 .  
## balance                     -6.487e-04  6.748e-04  -0.961 0.336393    
## shot_power                  -8.242e-04  6.693e-04  -1.231 0.218235    
## jumping                     -2.772e-04  4.665e-04  -0.594 0.552348    
## stamina                      4.075e-03  5.562e-04   7.326 2.47e-13 ***
## strength                     1.597e-04  6.561e-04   0.243 0.807728    
## long_shots                  -7.718e-04  7.095e-04  -1.088 0.276655    
## aggression                  -9.097e-04  5.067e-04  -1.795 0.072631 .  
## interceptions                7.142e-04  7.226e-04   0.988 0.322936    
## positioning                  5.412e-04  7.098e-04   0.763 0.445760    
## vision                       1.383e-03  6.627e-04   2.087 0.036870 *  
## penalties                    1.828e-05  6.238e-04   0.029 0.976629    
## composure                   -3.584e-05  7.386e-04  -0.049 0.961303    
## marking                      1.281e-03  5.884e-04   2.177 0.029521 *  
## standing_tackle             -2.176e-03  1.077e-03  -2.021 0.043266 *  
## sliding_tackle              -1.745e-03  1.003e-03  -1.739 0.082007 .  
## gk_diving                   -2.558e-03  1.385e-03  -1.847 0.064759 .  
## gk_handling                 -4.740e-03  1.394e-03  -3.400 0.000675 ***
## gk_kicking                  -1.530e-03  1.291e-03  -1.185 0.236041    
## gk_positioning              -1.812e-03  1.346e-03  -1.346 0.178178    
## gk_reflexes                 -1.730e-03  1.369e-03  -1.264 0.206102    
## release_clause               4.849e-01  8.764e-04 553.305  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5702 on 16585 degrees of freedom
## Multiple R-squared:  0.9901, Adjusted R-squared:  0.9901 
## F-statistic: 2.909e+04 on 57 and 16585 DF,  p-value: < 2.2e-16

Nhận xét

  • Dựa trên kết quả từ hàm summary, ta nhận thấy được độ chính xác của mô hình rất cao R-squared: 0.9901, Adjusted R-squared: 0.9901.
  • Với mức ý nghĩa alpha = 0.0001 và kết quả từ cột Pr(>|t|), các biến giải thích có ảnh hưởng đáng kể đến mô hình bao gồm: age, overall, potential, wage, international_reputation, positionalGoalKeeper, crossing, finishing, volleys, fk_accuracy, stamina, gk_kicking, release_clause. Tuy nhiên, cần lưu ý rằng các giá trị p-value này dựa trên các giả định của mô hình hồi quy tuyến tính.
  • Để đảm bảo tính chính xác và độ tin cậy của các hệ số ước lượng, nên áp dụng phương pháp bootstrap để kiểm định các hệ số thay vì chỉ dựa vào giả định ban đầu của mô hình.

1. Sự suy luận cho mô hình

  • Vì dữ liệu có kích thước rất lớn, thời gian thực hiện Boostrap sẽ rất lâu, nên ta thiết lập tham số parallel = multicore để thực hiện song song.
# Thực hiện phương pháp boostrap
fun_boot_md <- function(data, ind, formula, ...){
  data_new <- data[ind,]
  out_md <- lm(formula = formula, data = data_new, ...)
  return(setNames(out_md$coefficients, names(out_md$coefficients)))
}

set.seed(84)
boot_value_lm <- boot(data = simplified_data, statistic = fun_boot_md, R = 1000, formula = value ~ ., parallel = "multicore")

1.1 Kiểm định hệ số bằng phương pháp Bootstrap

pvals_boot_value <- sapply(1:ncol(boot_value_lm$t),function(x) {
  qt0 <- mean(boot_value_lm$t[, x] <= 0)
  if (qt0 < 0.5) {
    return(2*qt0)
  } else {
    return(2*(1- qt0))
  }
})
names(pvals_boot_value) <-  names(boot_value_lm$t0)
pvals_boot_value[pvals_boot_value < 0.0001]
##                      age                  overall                potential 
##                        0                        0                        0 
##                     wage international_reputation       positionGoalkeeper 
##                        0                        0                        0 
##                 crossing                finishing                  volleys 
##                        0                        0                        0 
##             long_passing                  stamina              gk_handling 
##                        0                        0                        0 
##           release_clause 
##                        0

Nhận xét

  • Với mức ý nghĩa alpha = 0.0001 các biến giải thích có ý nghĩa thống kê đối với mô hình bao gồm: age, overall, potential, international_reputation, positionGoalkeeper, crossing, finishing, short_passing, volleys, stamina, release_clause. Điều này cho thấy các biến này đều có ý nghĩa thống kê rất cao
  • Việc có nhiều biến có ý nghĩa thống kê cao thường cho thấy mô hình hồi quy được xây dựng là phù hợp và có khả năng giải thích tốt dữ liệu. Mặc dù các biến này đều có ý nghĩa thống kê, cần cân nhắc về ý nghĩa thực tiễn của chúng. Một số biến có thể không có ý nghĩa thực tế mạnh dù có ý nghĩa thống kê cao.

1.2 Đánh giá mô hình

set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_value_model <- train(value ~ ., data = simplified_data, method = "lm", trControl = ctrl)
print(cv_value_model)
## Linear Regression 
## 
## 16643 samples
##    52 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 14978, 14980, 14977, 14978, 14980, 14979, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.5730172  0.9901131  0.2563575
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Nhận xét:

  • Dựa vào kết quả từ đoạn code trên, với các giá trị RRMSE = 0.575626, R-squared = 0.573021, và MAE = 0.2563513, so sánh với kết quả của hàm summary.lm trong mục Mô hình hồi quy tuyến tính cho cột value, ta nhận thấy mô hình có khả năng giải thích rất cao khi sử dụng tất cả các biến độc lập. Điều này được thể hiện qua giá trị R-squared gần bằng 1, cho thấy mô hình phù hợp với dữ liệu và có độ chính xác cao trong việc dự đoán.
  • Tuy nhiên, mô hình hiện tại sử dụng rất nhiều biến giải thích, do đó, ta tìm cách đơn giản hóa mô hình

2. Lựa chọn mô hình

2.1 Hồi quy từng bước và cross-validation

predict.regsubsets <- function(object, newdata, id_model){
  form <- as.formula(object$call[[2]])
  x_mat <- model.matrix(form, newdata)
  coef_est <- coef(object, id = id_model)
  x_vars <- names(coef_est)
  res <- x_mat[, x_vars] %*% coef_est
  return(as.numeric(res))
}

# Chia bộ dữ liệu thành các folds
nrows <- nrow(simplified_data)
k <- 10

set.seed(21)
folds <-  sample(rep(1:k,length=nrows))

cv_error_rj <-  matrix(0,nrow=k,ncol=ncol(simplified_data))

for(r in 1:k){
  data_train_r  <-  simplified_data[folds != r,]
  data_test_r <-  simplified_data[folds == r,]
  steps_wise <- regsubsets(x=value ~ . ,data=data_train_r,  method="forward", really.big = TRUE, nvmax = ncol(simplified_data))
  
  for(j in 1:ncol(simplified_data)){
    pred_rj <-  predict(steps_wise, newdata = data_test_r,id_model = j)
    cv_error_rj[r,j]<-sqrt(mean((data_test_r$value - pred_rj)^2))
  }
}

cv_error <- colMeans(cv_error_rj)


ggplot(data= data.frame(x= c(1:ncol(simplified_data)),y=cv_error), mapping=aes(x=x,y=y)) +
  geom_point() +
  geom_line()+
  labs(x="Number of predictors",y="RMSE") +
  theme_bw()

Nhận xét:

  • Ta lựa chọn số biến cho mô hình tối ưu có số biến là 12
steps_wise <- regsubsets(x = value~ ., data=simplified_data, method="forward", really.big = TRUE, nvmax = ncol(simplified_data))

selected_steps_wise <- coef(steps_wise, id =  12)
selected_steps_wise
##              (Intercept)                      age                  overall 
##             -0.630166142             -0.018947474              0.031994567 
##                potential                     wage international_reputation 
##             -0.019646451              0.006282471              0.401685586 
##                 crossing                finishing                  volleys 
##             -0.002041508             -0.003374028              0.003115535 
##              fk_accuracy                  stamina           sliding_tackle 
##              0.001716851              0.004330145             -0.002808977 
##           release_clause 
##              0.485488417
# Tạo công thức cho mô hình hồi quy
steps_wise_formula <- as.formula(paste("value ~", paste(names(coef(steps_wise, id = 12))[-1], collapse = " + ")))
steps_wise_formula <- update(steps_wise_formula, . ~ . - positionForward - positionMidfielder)
steps_wise_formula <- update(steps_wise_formula, . ~ . + position)
steps_wise_formula
## value ~ age + overall + potential + wage + international_reputation + 
##     crossing + finishing + volleys + fk_accuracy + stamina + 
##     sliding_tackle + release_clause + position
# Sử dụng k-folds cross-validation để đánh giá mô hình
set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_steps_wise_value_model <- train(steps_wise_formula, data = simplified_data, method = "lm", trControl = ctrl)
print(cv_steps_wise_value_model)
## Linear Regression 
## 
## 16643 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 14978, 14980, 14977, 14978, 14980, 14979, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.5736017  0.9900932  0.2543194
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Nhận xét:

  • Sau khi giảm số biến giải thích xuống còn 12 biến, mô hình vẫn có khả năng giải thích rất cao với MSE = 0.5736017 và MAE = 0.2543194. Điều này cho thấy mô hình vẫn giữ được độ chính xác cao trong việc dự đoán giá trị của biến phụ thuộc, mặc dù số lượng biến giải thích đã giảm đáng kể. Việc giảm số biến giải thích giúp đơn giản hóa mô hình, làm cho nó dễ hiểu hơn và giảm nguy cơ overfitting (quá khớp) khi áp dụng trên dữ liệu mới.

2.2 Hồi quy lasso

Lasso (Least Absolute Shrinkage and selection Operator) nổi tiếng với khả năng lựa chọn biến bằng cách co một số hệ số về 0. Điều này rất hữu ích khi xử lý dữ liệu có nhiều biến dự đoán. Tuy nhiên, Lasso vẫn giả định mối quan hệ tuyến tính giữa các biến dự đoán còn lại và biến kết quả.

x_data <- model.matrix(value ~ ., data = simplified_data)[,-1]
y_data <- simplified_data$value


# Tìm hệ só lambda tối ưu bằng phương pháp cross-validation
set.seed(24)
out_cv_lasso_value <- cv.glmnet(x = x_data, y = y_data, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_value <- out_cv_lasso_value$lambda.min
cat("Hệ số lambda tối ưu là: ", lambda_lasso_value)
## Hệ số lambda tối ưu là:  0.02579507
# Từ giá trị lambda tối ưu tìm được, fit mô hình trên toàn bộ tập data
out_lasso_md <- glmnet(x = x_data, y = y_data, alpha = 1, lambda = lambda_lasso_value, family = "gaussian")

# Các hệ số trong mô hình
coeff_lasso_md <- predict(out_lasso_md, type = "coefficients")
nonzero_coeff_lasso_md <- setNames(coeff_lasso_md@x, rownames(coeff_lasso_md)[coeff_lasso_md@i + 1])
nonzero_coeff_lasso_md
##              (Intercept)                  overall                     wage 
##            -1.1648293821             0.0134804901             0.0055903711 
## international_reputation                  volleys                  stamina 
##             0.3351279874             0.0006692361             0.0004724311 
##           release_clause 
##             0.4863211056

Nhận xét

  • Các biến dự đoán có ý nghĩa: Các hệ số của các biến dự đoán sau không bằng không, cho thấy chúng có ý nghĩa trong việc dự đoán biến phản hồi: overall: 0.0131, wage: 5.6743, international_reputation: 0.3363, volleys: 0.0007, reactions: 0.0002, stamina: 0.0005, release_clause: 0.4862
  • Diễn giải:
    • wage có hệ số dương lớn nhất, cho thấy nó có tác động mạnh mẽ đến biến phản hồi.
    • release_clause cũng có hệ số dương đáng kể, cho thấy nó là một biến dự đoán quan trọng.
    • Các hệ số không bằng không của overall, international_reputation, volleys, reactions, và stamina chỉ ra rằng các biến này cũng đóng góp vào mô hình, mặc dù ít hơn.

Đánh giá mô hình

set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))

cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)

for (i in 1:k) {
  train_indices <- which(folds != i)
  test_indices <- which(folds == i)
  
  train_data <- simplified_data[train_indices, ]
  test_data <- simplified_data[test_indices, ]
  
  x_train <- model.matrix(value ~ ., data = train_data)[, -1]
  y_train <- train_data$value
  x_test <- model.matrix(value ~ ., data = test_data)[, -1]
  y_test <- test_data$value
  
  out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_value, family = "gaussian")
  pred <- predict(out_lasso, s = lambda_lasso_value, newx = x_test)
  
  cv_rmse[i] <- rmse(y_test, pred)
  cv_rsquare[i] <- cor(y_test, pred)^2
  cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}

avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)

cat("Average RMSE: ", avg_rmse)
## Average RMSE:  0.5781515
cat("\nAverage R-squared: ", avg_rsquare)
## 
## Average R-squared:  0.9897874
cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)
## 
## Average R-squared adjusted:  0.9897874

3 Chuẩn đoán mô hình

3.1 Kiểm tra tính tuyến tính và Tính đồng nhất phương sai

pred_lasso_md <- predict(out_lasso_md, newx = x_data)
residual_lasso_md <- pred_lasso_md - simplified_data$value

ggplot(data = simplified_data, mapping = aes(x = pred_lasso_md, y = residual_lasso_md)) +
 geom_point() +
 geom_smooth(method = "loess", se = FALSE) +
 geom_hline(yintercept = 0, linetype = "dashed") +
 labs(x = "Fittted values", y = "Residuals") +
 theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét

  • Đường trung bình phần dư (đường màu xanh) khá tương đồng với đường thẳng. Điều này có thể là dấu hiệu của mối quan hệ tuyến tính giữa value và một hoặc nhiều biến hồi quy.
  • Các cầu thủ có giá trị thấp chiếm số lượng lớn, gây ra sự tập trung dày đặc trong đồ thị phần dư ở khu vực fitted values nhỏ. Ngược lại, số lượng cầu thủ có giá trị cao ít hơn, dẫn đến phần dư phân tán nhiều hơn ở các giá trị fitted values lớn. Điều này cho thấy phương sai của phần dư không đồng đều, vi phạm giả định về homoscedasticity (phương sai phần dư đồng đều) trong hồi quy tuyến tính.

3.2 Kiểm tra tính tuyến tính từng phần

fitted_lasso_md <- predict(out_lasso_md, newx = x_data)
resid_lasso_md <- y_data - fitted_lasso_md
selected_var_lasso_md <- names(nonzero_coeff_lasso_md[-1, drop = FALSE])


for(col_name in selected_var_lasso_md){
    terms <- x_data[, col_name] * nonzero_coeff_lasso_md[col_name]
    p <- ggplot(x_data, mapping = aes(x_data[, col_name], resid_lasso_md + terms)) +
    geom_point() +
    labs(x = col_name, y = "Partial Residuals") +
    geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
    geom_line(aes(x = x_data[, col_name], y = terms), color = "blue")
    theme_bw()
    print(p)
}
## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

Nhận xét

overall:
- Mối quan hệ phi tuyến tính: Đường smooth có xu hướng tăng lên khi giá trị “overall” tăng lên. Điều này cho thấy mối quan hệ giữa biến “overall” và biến kết quả có thể không tuyến tính - Heteroscedasticity (phương sai không đồng nhất): Sự phân tán của các điểm quanh đường smooth không đồng đều. Khi giá trị “overall” tăng lên (khoảng 75 trở lên), các điểm trở nên phân tán rộng hơn, đặc biệt là ở các giá trị trên 80.

international_reputation:
- Mối quan hệ tuyến tính yếu: “international_reputation” trong mô hình Lasso. Đường màu xanh dương có độ dốc rất nhỏ, gần như bằng phẳng, cho thấy mối quan hệ tuyến tính giữa biến này và biến kết quả là rất yếu.
- Heteroscedasticity: Sự phân tán của các điểm có vẻ không đồng đều giữa các mức độ. Có vẻ như ở mức độ 3, các điểm phân tán rộng hơn so với các mức độ khác.

Voleys:
- Mối quan hệ phi tuyến tính nhẹ: Đường smooth cho thấy một mối quan hệ phi tuyến tính nhẹ giữa “volleys” và biến kết quả. Ở các giá trị “volleys” thấp và trung bình, đường smooth khá bằng phẳng, nhưng có xu hướng hơi tăng ở các giá trị “volleys” cao. - Phân bố dữ liệu tập trung: Phần lớn dữ liệu tập trung ở các giá trị “volleys” từ khoảng 0 đến 80. release_clause:
- Mối quan hệ tuyến tính mạnh mẽ: Đường smooth gần như trùng khớp với đường thẳng màu xanh dương, cho thấy một mối quan hệ tuyến tính rất mạnh mẽ giữa “release_clause” và biến kết quả. - Sự phân tán tương đối đồng đều: Các điểm dữ liệu phân tán khá đều quanh đường thẳng, không có dấu hiệu rõ ràng của heteroscedasticity (phương sai không đồng nhất).

4. Mở rộng mô hình

4.1 Hồi quy đa thức

  • Dựa vào các biểu đồ kiểm tra tính tuyến tính từng phần, ta thấy rằng các biến overall, wage, volleysrelease_clause có mối quan hệ phi tuyến tính với biến phụ thuộc. Do đó, ta sẽ xây dựng mô hình hồi quy đa thức cho các biến này.
poly_value_md <- lm(value ~  release_clause + poly(overall, 3) + poly(wage, 3) + international_reputation + poly(volleys,2) + stamina + release_clause*wage, data = simplified_data)

summary(poly_value_md)
## 
## Call:
## lm(formula = value ~ release_clause + poly(overall, 3) + poly(wage, 
##     3) + international_reputation + poly(volleys, 2) + stamina + 
##     release_clause * wage, data = simplified_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7527 -0.1113 -0.0088  0.1034 10.8366 
## 
## Coefficients: (1 not defined because of singularities)
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.335e-01  3.072e-02   7.600 3.12e-14 ***
## release_clause            4.051e-01  1.509e-03 268.366  < 2e-16 ***
## poly(overall, 3)1         8.089e+01  1.462e+00  55.331  < 2e-16 ***
## poly(overall, 3)2         7.567e+01  1.323e+00  57.214  < 2e-16 ***
## poly(overall, 3)3         3.804e+01  8.405e-01  45.255  < 2e-16 ***
## poly(wage, 3)1           -5.363e+01  1.909e+00 -28.101  < 2e-16 ***
## poly(wage, 3)2           -5.573e+01  1.532e+00 -36.378  < 2e-16 ***
## poly(wage, 3)3           -6.201e+00  5.747e-01 -10.790  < 2e-16 ***
## international_reputation -1.051e-01  1.626e-02  -6.463 1.06e-10 ***
## poly(volleys, 2)1         6.634e+00  6.413e-01  10.345  < 2e-16 ***
## poly(volleys, 2)2         1.264e+01  6.449e-01  19.602  < 2e-16 ***
## stamina                   5.924e-03  3.522e-04  16.821  < 2e-16 ***
## wage                             NA         NA      NA       NA    
## release_clause:wage       3.704e-04  9.624e-06  38.483  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5115 on 16630 degrees of freedom
## Multiple R-squared:  0.992,  Adjusted R-squared:  0.992 
## F-statistic: 1.721e+05 on 12 and 16630 DF,  p-value: < 2.2e-16

Đánh giá mô hình bằng phương pháp cross-validation

set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))

cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)

for (i in 1:k) {
  train_indices <- which(folds != i)
  test_indices <- which(folds == i)
  
  train_data <- simplified_data[train_indices, ]
  test_data <- simplified_data[test_indices, ]
  y_test <- test_data$value
  
  poly_value_md <- lm(value ~  release_clause + poly(overall, 3) + poly(wage, 3) + international_reputation + poly(volleys,2) + stamina , data = train_data)
  pred <- predict(poly_value_md, newdata = test_data)
  

  cv_rmse[i] <- rmse(y_test, pred)
  cv_rsquare[i] <- cor(y_test, pred)^2
  cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}


avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)

cat("Average RMSE: ", avg_rmse)
## Average RMSE:  0.5401948
cat("\nAverage R-squared: ", avg_rsquare)
## 
## Average R-squared:  0.9911411
cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)
## 
## Average R-squared adjusted:  0.9911411

4.2 Mô hình hồi quy tổng quát

gam_value_md <- gam(value ~  release_clause + s(overall) + s(wage) + international_reputation + s(volleys) + stamina + release_clause*wage, data = simplified_data)
summary(gam_value_md)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## value ~ release_clause + s(overall) + s(wage) + international_reputation + 
##     s(volleys) + stamina + release_clause * wage
## 
## Parametric coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.014e+00  5.044e-02  20.101  < 2e-16 ***
## release_clause            4.083e-01  1.563e-03 261.254  < 2e-16 ***
## international_reputation -1.302e-01  1.603e-02  -8.124 4.83e-16 ***
## stamina                   4.692e-03  3.679e-04  12.754  < 2e-16 ***
## wage                     -7.031e-02  4.268e-03 -16.471  < 2e-16 ***
## release_clause:wage       3.197e-04  1.134e-05  28.196  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F p-value    
## s(overall) 8.932  8.997 474.8  <2e-16 ***
## s(wage)    7.932  8.140 223.5  <2e-16 ***
## s(volleys) 8.978  9.000 115.6  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 32/33
## R-sq.(adj) =  0.993   Deviance explained = 99.3%
## GCV =  0.241  Scale est. = 0.24054   n = 16643
predictions <- predict(gam_value_md, newdata = simplified_data)
mse <- mean((simplified_data$value - predictions)^2)
rmse <- sqrt(mse)
print(rmse)
## [1] 0.4899834

Đánh giá mô hình bằng phương pháp cross-validation

set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))

cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)

for (i in 1:k) {
  train_indices <- which(folds != i)
  test_indices <- which(folds == i)
  
  train_data <- simplified_data[train_indices, ]
  test_data <- simplified_data[test_indices, ]
  y_test <- test_data$value
  
  gam_value_md <- gam(value ~  release_clause + s(overall) + s(wage) + international_reputation + s(volleys) + stamina + release_clause*wage, data = train_data)
  pred <- predict(gam_value_md, newdata = test_data)
  

  cv_rmse[i] <- rmse(y_test, pred)
  cv_rsquare[i] <- cor(y_test, pred)^2
  cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}


avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)

cat("Average RMSE: ", avg_rmse)
## Average RMSE:  0.5166711
cat("\nAverage R-squared: ", avg_rsquare)
## 
## Average R-squared:  0.9917792
cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)
## 
## Average R-squared adjusted:  0.9917792

5. Kết luận

  • Sau khi xây dựng mô hình hồi quy tuyến tính, hồi quy lasso, hồi quy đa thức và hồi quy GAM, ta nhận thấy rằng mô hồi quy GAM có khả năng giải thích tốt nhất với Average RMSE: 0.5166711 Average R-squared: 0.9917792, Average R-squared adjusted: 0.9917792.
  • Tổng số biến giải thích sử dụng cho mô hình này là 6 biến giải thích, bao gồm: release_clause, overall, wage, international_reputation, volleysstamina.

IV. Linear Regresion cho wage

wage_lm <- lm(wage ~ ., data = simplified_data)

summary(wage_lm)
## 
## Call:
## lm(formula = wage ~ ., data = simplified_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -154.819   -2.072   -0.239    1.581  195.650 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  7.183e+02  1.566e+02   4.587 4.53e-06 ***
## age                          2.732e-01  4.095e-02   6.671 2.62e-11 ***
## overall                     -3.915e-02  4.690e-02  -0.835 0.403933    
## potential                    2.165e-02  3.402e-02   0.637 0.524399    
## value                        2.232e+00  1.443e-01  15.471  < 2e-16 ***
## preferred_footRight          3.728e-02  2.086e-01   0.179 0.858152    
## international_reputation     9.678e+00  2.995e-01  32.317  < 2e-16 ***
## weak_foot                   -1.681e-01  1.379e-01  -1.219 0.222699    
## skill_moves                 -9.454e-01  2.256e-01  -4.190 2.80e-05 ***
## attacking_work_rateLow      -4.702e-01  4.538e-01  -1.036 0.300182    
## attacking_work_rateMedium    5.003e-03  2.178e-01   0.023 0.981674    
## deffensive_work_rate Low    -1.197e+00  3.941e-01  -3.037 0.002394 ** 
## deffensive_work_rate Medium -9.628e-01  2.449e-01  -3.932 8.47e-05 ***
## body_typeNormal             -3.246e-01  1.942e-01  -1.671 0.094706 .  
## body_typeStocky             -8.245e-01  3.971e-01  -2.076 0.037880 *  
## positionForward              1.254e+00  4.772e-01   2.627 0.008614 ** 
## positionGoalkeeper           4.857e+00  1.799e+00   2.700 0.006937 ** 
## positionMidfielder          -2.817e-01  3.358e-01  -0.839 0.401551    
## jersey_number                1.412e-02  5.465e-03   2.585 0.009758 ** 
## joined                      -2.529e-01  4.155e-02  -6.085 1.19e-09 ***
## contract_valid_until        -1.167e-01  6.944e-02  -1.680 0.092934 .  
## height                       1.357e-01  6.489e-02   2.091 0.036537 *  
## weight                       1.939e-02  9.822e-03   1.974 0.048385 *  
## crossing                     3.842e-02  1.186e-02   3.241 0.001194 ** 
## finishing                   -2.069e-02  1.439e-02  -1.437 0.150666    
## heading_accuracy             1.683e-02  1.267e-02   1.328 0.184249    
## short_passing               -1.246e-02  2.030e-02  -0.614 0.539470    
## volleys                     -9.718e-03  1.229e-02  -0.791 0.429220    
## dribbling                    1.899e-02  1.778e-02   1.068 0.285683    
## curve                        4.183e-03  1.213e-02   0.345 0.730257    
## fk_accuracy                 -3.701e-02  1.077e-02  -3.436 0.000591 ***
## long_passing                -2.351e-02  1.494e-02  -1.573 0.115647    
## ball_control                 4.969e-02  2.216e-02   2.242 0.024945 *  
## acceleration                -1.537e-02  1.664e-02  -0.924 0.355628    
## sprint_speed                 2.035e-02  1.544e-02   1.318 0.187380    
## agility                     -1.337e-05  1.235e-02  -0.001 0.999136    
## reactions                   -3.647e-02  1.883e-02  -1.937 0.052771 .  
## balance                      3.492e-02  1.262e-02   2.766 0.005678 ** 
## shot_power                   2.591e-02  1.252e-02   2.069 0.038567 *  
## jumping                      1.110e-03  8.728e-03   0.127 0.898806    
## stamina                     -4.060e-02  1.042e-02  -3.897 9.78e-05 ***
## strength                    -2.356e-02  1.227e-02  -1.920 0.054919 .  
## long_shots                   1.825e-03  1.327e-02   0.137 0.890667    
## aggression                  -6.621e-04  9.482e-03  -0.070 0.944329    
## interceptions                1.345e-02  1.352e-02   0.995 0.319690    
## positioning                  2.545e-03  1.328e-02   0.192 0.848007    
## vision                       4.666e-03  1.240e-02   0.376 0.706754    
## penalties                    3.875e-02  1.167e-02   3.321 0.000899 ***
## composure                   -1.655e-02  1.382e-02  -1.198 0.231026    
## marking                     -1.717e-02  1.101e-02  -1.560 0.118850    
## standing_tackle             -8.809e-03  2.015e-02  -0.437 0.661933    
## sliding_tackle               7.356e-02  1.876e-02   3.920 8.88e-05 ***
## gk_diving                    5.206e-04  2.591e-02   0.020 0.983971    
## gk_handling                  3.398e-02  2.609e-02   1.303 0.192762    
## gk_kicking                  -1.186e-02  2.415e-02  -0.491 0.623374    
## gk_positioning              -5.471e-02  2.518e-02  -2.173 0.029793 *  
## gk_reflexes                 -6.646e-03  2.561e-02  -0.260 0.795239    
## release_clause               3.560e-01  7.228e-02   4.925 8.51e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.67 on 16585 degrees of freedom
## Multiple R-squared:  0.7711, Adjusted R-squared:  0.7703 
## F-statistic: 980.4 on 57 and 16585 DF,  p-value: < 2.2e-16

1. Lựa chọn mô hình

1.1 Hồi quy từng bước và cross-validation

predict.regsubsets <- function(object, newdata, id_model){
  form <- as.formula(object$call[[2]])
  x_mat <- model.matrix(form, newdata)
  coef_est <- coef(object, id = id_model)
  x_vars <- names(coef_est)
  res <- x_mat[, x_vars] %*% coef_est
  return(as.numeric(res))
}

# Chia bộ dữ liệu thành các folds
nrows <- nrow(simplified_data)
k <- 10

set.seed(23)
folds <-  sample(rep(1:k,length=nrows))

cv_error_rj <-  matrix(0,nrow=k,ncol=ncol(simplified_data))

for(r in 1:k){
  data_train_r  <-  simplified_data[folds != r,]
  data_test_r <-  simplified_data[folds == r,]
  steps_wise <- regsubsets(x=wage ~ . ,data=data_train_r,  method="forward", really.big = TRUE, nvmax = ncol(simplified_data))
  
  for(j in 1:ncol(simplified_data)){
    pred_rj <-  predict(steps_wise, newdata = data_test_r,id_model = j)
    cv_error_rj[r,j]<-sqrt(mean((data_test_r$value - pred_rj)^2))
  }
}

cv_error <- colMeans(cv_error_rj)

ggplot(data= data.frame(x= c(1:ncol(simplified_data)),y=cv_error), mapping = aes(x = x,y=y)) +
  geom_point() +
  geom_line()+
  labs(x="Number of predictors",y="RMSE") +
  theme_bw()

out_subset_wage <- regsubsets(x = wage ~ ., data = simplified_data, method = "forward", nvmax = ncol(simplified_data))
coef(out_subset_wage, id = which.min(cv_error))
## (Intercept)       value 
##    1.428039    3.352893

Đánh giá mô hình

# Tạo công thức cho mô hình hồi quy
subset_wage_formula <- as.formula(paste("wage ~", paste(names(coef(out_subset_wage, id = which.min(cv_error)))[-1], collapse = " + ")))
#subset_wage_formula <- update(subset_wage_formula, . ~ . - positionMidfielder)
#subset_wage_formula <- update(subset_wage_formula, . ~ . + position)
subset_wage_formula
## wage ~ value
# Sử dụng k-folds cross-validation để đánh giá mô hình
set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_subset_wage <- train(subset_wage_formula, data = simplified_data, method = "lm", trControl = ctrl)
print(cv_subset_wage)
## Linear Regression 
## 
## 16643 samples
##     1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 14978, 14980, 14978, 14979, 14980, 14979, ... 
## Resampling results:
## 
##   RMSE      Rsquared  MAE     
##   11.27782  0.736994  4.767674
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

1.2 Hồi quy lasso

Lasso (Least Absolute Shrinkage and selection Operator) nổi tiếng với khả năng lựa chọn biến bằng cách co một số hệ số về 0. Điều này rất hữu ích khi xử lý dữ liệu có nhiều biến dự đoán. Tuy nhiên, Lasso vẫn giả định mối quan hệ tuyến tính giữa các biến dự đoán còn lại và biến kết quả.

x_data_wage <- model.matrix(wage ~ ., data = simplified_data)[,-1]
y_data_wage <- simplified_data$wage


# Tìm hệ só lambda tối ưu bằng phương pháp cross-validation
set.seed(21)
out_cv_lasso_wage <- cv.glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_wage <- out_cv_lasso_wage$lambda.min
cat("Hệ số lambda tối ưu là: ", lambda_lasso_wage)
## Hệ số lambda tối ưu là:  0.01353117
# Từ giá trị lambda tối ưu tìm được, fit mô hình trên toàn bộ tập data
out_lasso_wage <- glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")

# Các hệ số trong mô hình
coeff_lasso_wage <- predict(out_lasso_wage, type = "coefficients")
nonzero_coeff_lasso_wage <- setNames(coeff_lasso_wage@x, rownames(coeff_lasso_wage)[coeff_lasso_wage@i + 1])
sort(nonzero_coeff_lasso_wage, decreasing = T)
##                 (Intercept)    international_reputation 
##                6.958199e+02                9.720065e+00 
##          positionGoalkeeper                       value 
##                2.720660e+00                2.321381e+00 
##             positionForward              release_clause 
##                1.081515e+00                3.056566e-01 
##                         age                      height 
##                2.357464e-01                1.136092e-01 
##              sliding_tackle                    crossing 
##                6.288031e-02                3.665628e-02 
##                ball_control                   penalties 
##                3.308576e-02                3.252654e-02 
##                     balance                  shot_power 
##                2.718679e-02                2.111946e-02 
##                      weight                   dribbling 
##                1.539182e-02                1.525324e-02 
##                 gk_handling               jersey_number 
##                1.392103e-02                1.384133e-02 
##            heading_accuracy               interceptions 
##                1.033243e-02                5.719712e-03 
##                sprint_speed                       curve 
##                3.291576e-03                1.283051e-04 
##                     jumping                acceleration 
##                9.939096e-05               -3.373238e-04 
##               short_passing                     volleys 
##               -7.539697e-04               -3.246271e-03 
##                   finishing                     marking 
##               -1.013665e-02               -1.106648e-02 
##                   composure                     overall 
##               -1.118500e-02               -1.766929e-02 
##                    strength                long_passing 
##               -1.879400e-02               -1.953088e-02 
##              gk_positioning                   reactions 
##               -2.442599e-02               -2.965283e-02 
##                 fk_accuracy                     stamina 
##               -3.175789e-02               -3.775269e-02 
##        contract_valid_until                   weak_foot 
##               -1.083770e-01               -1.463178e-01 
##                      joined             body_typeNormal 
##               -2.483883e-01               -2.715708e-01 
##          positionMidfielder      attacking_work_rateLow 
##               -4.079299e-01               -4.415127e-01 
##             body_typeStocky deffensive_work_rate Medium 
##               -7.319435e-01               -8.779270e-01 
##                 skill_moves    deffensive_work_rate Low 
##               -8.874190e-01               -1.075323e+00

Nhận xét

  • Từ hệ số của các biến trong mô hình hồi quy Lasso, ta thấy chỉ có hệ số của các biến giải thích như overall, international_reputation, volleys, reactions, stamina, release_clause khác 0, cho thấy các biến này có ảnh hưởng đến biến kết quả wage.
  • Từ hệ số của các biến trong mô hình hồi quy Lasso, ta các biến giải thích sau international_reputation, position, value, release_clause có hệ số lớn nhât

Đánh giá mô hình

set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))

cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)

for (i in 1:k) {
  train_indices <- which(folds != i)
  test_indices <- which(folds == i)
  
  train_data <- simplified_data[train_indices, ]
  test_data <- simplified_data[test_indices, ]
  
  x_train <- model.matrix(wage ~ ., data = train_data)[, -1]
  y_train <- train_data$wage
  x_test <- model.matrix(wage ~ ., data = test_data)[, -1]
  y_test <- test_data$wage
  
  out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
  pred <- predict(out_lasso, s = lambda_lasso_wage, newx = x_test)
  
  cv_rmse[i] <- rmse(y_test, pred)
  cv_rsquare[i] <- cor(y_test, pred)^2
  cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}

avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)

cat("Average RMSE: ", avg_rmse)
## Average RMSE:  10.664
cat("\nAverage R-squared: ", avg_rsquare)
## 
## Average R-squared:  0.7649655
cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)
## 
## Average R-squared adjusted:  0.7649655
  • Từ kết quả của mô hình lasso, ta đánh giá mô hình chỉ dựa trên các biến giải thích có hệ số lớn như: international_reputation, position, value, deffensive_work_rate, skill_moves
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))

cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)

for (i in 1:k) {
  train_indices <- which(folds != i)
  test_indices <- which(folds == i)
  
  train_data <- simplified_data[train_indices, ]
  test_data <- simplified_data[test_indices, ]
  
  x_train <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = train_data)[, -1]
  y_train <- train_data$wage
  x_test <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = test_data)[, -1]
  y_test <- test_data$wage
  
  out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
  pred <- predict(out_lasso, s = lambda_lasso_wage, newx = x_test)
  
  cv_rmse[i] <- rmse(y_test, pred)
  cv_rsquare[i] <- cor(y_test, pred)^2
  cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}

avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)

cat("Average RMSE: ", avg_rmse)
## Average RMSE:  10.74997
cat("\nAverage R-squared: ", avg_rsquare)
## 
## Average R-squared:  0.7610392
cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)
## 
## Average R-squared adjusted:  0.7610392

Nhận xét

  • Mô hình Lasso sử dụng các biến giải thích với hệ số khác 0 cho kết quả: Average RMSE: 10.66493, Average R-squared: 0.7649227, Average R-squared adjusted: 0.7649227. Điều này cho thấy mô hình không có khả năng giải thích tốt
  • Sau khi chọn ra các biến giải thích như: international_reputation, position, value, deffensive_work_rate, skill_moves, mô hình có khả năng giải thích tốt hơn với Average RMSE: 10.74997, Average R-squared: 0.7610392, Average R-squared adjusted: 0.7610392. Ta thấy cả 2 mô hình không có sự khác biệt lớn về độ chính xác giữa chúng. Do đó, ta sẽ chọn mô hình Lasso sử dụng ít biến giải thích hơn để giảm độ phức tạp của mô hình.
x_data_wage <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = simplified_data)[,-1]
y_data_wage <- simplified_data$wage

out_cv_lasso_wage <- cv.glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_wage <- out_cv_lasso_wage$lambda.min
out_lasso_wage <- glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
coeff_lasso_wage <- predict(out_lasso_wage, type = "coefficients")
nonzero_coeff_lasso_wage <- setNames(coeff_lasso_wage@x, rownames(coeff_lasso_wage)[coeff_lasso_wage@i + 1])
nonzero_coeff_lasso_wage
##                 (Intercept)    international_reputation 
##                  -6.9870834                  10.8537193 
##             positionForward          positionGoalkeeper 
##                  -0.4267046                  -1.1183289 
##          positionMidfielder                       value 
##                  -1.1458339                   2.8598020 
##    deffensive_work_rate Low deffensive_work_rate Medium 
##                  -1.6097663                  -1.2082737 
##                 skill_moves 
##                  -0.3456851

2. Chuẩn đoán mô hình

2.1 Kiểm tra tính tuyến tính và Tính đồng nhất phương sai

pred_lasso_wage <- predict(out_lasso_wage, newx = x_data_wage)

residual_lasso_wage <- pred_lasso_wage - simplified_data$wage

ggplot(data = simplified_data, mapping = aes(x = pred_lasso_wage, y = residual_lasso_wage)) +
 geom_point() +
 geom_smooth(method = "loess", se = FALSE) +
 geom_hline(yintercept = 0, linetype = "dashed") +
 labs(x = "Fittted values", y = "Residuals") +
 theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét

  • Đường trung bình phần dư (đường màu xanh) có xu hướng cong xuống khi giá trị biến phụ thuộc tăng lên. Điều này có thể là dấu hiệu của mối quan hệ phi tuyến tuyến tính giữa wage và các biến giải thích
  • Các cầu thủ có giá trị thấp chiếm số lượng lớn, gây ra sự tập trung dày đặc trong đồ thị phần dư ở khu vực fitted values nhỏ. Ngược lại, số lượng cầu thủ có giá trị cao ít hơn, dẫn đến phần dư phân tán nhiều hơn ở các giá trị fitted values lớn. Ở các cầu thủ có giá trị lớn, các điểm phân tán xa đường ngang residuals 0 hơn. Vì vậy, mô hình không đáp ứng được giả định về tính đồng nhất phương sai.

2.2 Kiểm tra tính tuyến tính từng phần

fitted_lasso_wage <- predict(out_lasso_wage, newx = x_data_wage)
resid_lasso_wage <- simplified_data$wage - fitted_lasso_wage

for(col_name in colnames(x_data_wage)){
    terms <- x_data_wage[, col_name] * nonzero_coeff_lasso_wage[col_name]
    p <- ggplot(x_data_wage, mapping = aes(x_data_wage[, col_name], resid_lasso_wage + terms)) +
    geom_point() +
    labs(x = col_name, y = "Partial Residuals") +
    geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
    geom_line(aes(x = x_data_wage[, col_name], y = terms), color = "blue")
    theme_bw()
    print(p)
}
## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

Nhận xét

international_reputation:
- Mối quan hệ tuyến tính mạnh mẽ: Đường smooth gần như trùng khớp với đường thẳng màu xanh dương, cho thấy một mối quan hệ tuyến tính rất mạnh mẽ giữa “international_reputation” và biến kết quả.

positionForward, positionGoalkeeper, positionMidfieler, defensive_work_rate, skill_moves:
- Các phần dư có sự phân tán lớn, đặc biệt ở cả hai đầu. Điều này có thể cho thấy một mức độ biến thiên đáng kể mà biến positionForward, positionaMildfielder không thể giải thích.

3. Mở rộng mô hình

3.1 Mô hình hồi quy đa thức

poly_wage_md <- lm(wage ~  international_reputation + poly(value, 3) +  poly(skill_moves, 3)  , data = simplified_data)

summary(poly_wage_md)
## 
## Call:
## lm(formula = wage ~ international_reputation + poly(value, 3) + 
##     poly(skill_moves, 3), data = simplified_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -163.952   -1.906   -0.955    0.844  191.303 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                -3.4788     0.3234 -10.756  < 2e-16 ***
## international_reputation   11.7479     0.2804  41.897  < 2e-16 ***
## poly(value, 3)1          2122.5095    15.2559 139.128  < 2e-16 ***
## poly(value, 3)2           145.6990    11.1851  13.026  < 2e-16 ***
## poly(value, 3)3           -35.8573    10.9297  -3.281  0.00104 ** 
## poly(skill_moves, 3)1     -19.5235    11.8291  -1.650  0.09887 .  
## poly(skill_moves, 3)2    -107.2838    11.3589  -9.445  < 2e-16 ***
## poly(skill_moves, 3)3     -83.2415    10.8345  -7.683 1.64e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.71 on 16635 degrees of freedom
## Multiple R-squared:  0.7689, Adjusted R-squared:  0.7688 
## F-statistic:  7905 on 7 and 16635 DF,  p-value: < 2.2e-16

3.2 Mô hình tuyến tính tổng quát

gam_wage_md <- gam(wage ~  international_reputation  + s(value) + value * international_reputation + skill_moves, data = simplified_data)
summary(gam_wage_md)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## wage ~ international_reputation + s(value) + value * international_reputation + 
##     skill_moves
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    14.09120    1.10488  12.754  < 2e-16 ***
## international_reputation        3.35370    0.36048   9.303  < 2e-16 ***
## value                          -4.59308    0.41214 -11.145  < 2e-16 ***
## skill_moves                    -0.35949    0.11642  -3.088  0.00202 ** 
## international_reputation:value  0.91223    0.02533  36.014  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##            edf Ref.df     F p-value    
## s(value) 8.159  8.175 164.2  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 13/14
## R-sq.(adj) =  0.789   Deviance explained =   79%
## GCV = 104.46  Scale est. = 104.38    n = 16643
predictions <- predict(gam_wage_md, newdata = simplified_data)
mse <- mean((simplified_data$wage - predictions)^2)
rmse <- sqrt(mse)
print(rmse)
## [1] 10.21275

4. Kết luận

  • So với mô hình ban đầu, mô hình cuối cùng chỉ sử dụng các biến international_reputation, value, skill_moves có kết quả tương đồng, nhưng mô hình cuối cùng đơn giản hơn rất nhiều
  • Kết quả từ mô hình với R-squared-adjusted = 0.789, và RMSE = 10.21275 cho thấy mô hình hồi quy cho biến giải thích wage không có khả năng giải thích cho bộ dữ liệu